##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineTpm.FilterLowExpression-MergeMutiSample.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/showGene"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"
    gene <- "CDX2"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
gtf_file <- opt$gtf_file
gene <- opt$gene

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_expression <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))

colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

igc_class <- c("IM + IGC + DGC" , "IM + IGC")
dgc_class <- c("IM + IGC + DGC" , "IM + DGC")

im_igc_sample <- unique(info[info$Type %in% igc_class , "ID"])
im_dgc_sample <- unique(info[info$Type %in% dgc_class , "ID"])

##########################################################################################
plotTpm <- function(tmp_dat = tmp_dat , out_name = out_name , title = title){

    my_comparisons_1 <- list( 
        c(1, 2), c(1, 3) , 
        c(2, 3) )

    plot <- ggplot( tmp_dat , aes( x = Class , y = Tpm , color = Class ) ) +
        geom_line( aes( group = Sample ) , size = 0.4 , color = "gray" ) +  ## 配对样本加线
        geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        facet_grid(.~Type,space='free_x',scales='free_x') +
        scale_fill_npg()+
        scale_color_npg()+
        xlab(NULL) +
        ylab("TPM")+
        theme_bw() +
        ggtitle(title) +
        stat_compare_means(comparisons = my_comparisons_1) +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 8,color="black",face='bold'),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5)) 
    ggsave(file=out_name,plot=plot,width=5,height=5)

}

describeGene <- function(dat_expression = dat_expression , image_path = image_path ){

    tmp_exp <- subset( dat_expression , Hugo_Symbol == gene)

    sample <- sapply( strsplit(colnames(tmp_exp)[3:(ncol(tmp_exp))] , "_") , "[" , 1)
    class <- sapply( strsplit(colnames(tmp_exp)[3:(ncol(tmp_exp))] , "_") , "[" , 2)
    tpm <- as.numeric(tmp_exp[3:(ncol(tmp_exp))])
    tmp_dat <- data.frame( Sample = sample , Class = class , Tpm = tpm )

    ## Normal->IM->IGC看基因的表达变化情况
    tmp_dat_igc <- subset( tmp_dat , Sample %in% im_igc_sample )
    tmp_dat_igc <- subset( tmp_dat_igc , Class %in% c( "Normal" , "IM" , "IGC" ))
    tmp_dat_igc$Class <- factor( tmp_dat_igc$Class , levels = c( "Normal" , "IM" , "IGC" ) , order = T )
    tmp_dat_igc$Type <- "IM + IGC"

    ## Normal->IM->DGC看基因的表达变化情况
    tmp_dat_dgc <- subset( tmp_dat , Sample %in% im_dgc_sample )
    tmp_dat_dgc <- subset( tmp_dat_dgc , Class %in% c( "Normal" , "IM" , "DGC" ))
    tmp_dat_dgc$Class <- factor( tmp_dat_dgc$Class , levels = c( "Normal" , "IM" , "DGC" ) , order = T )
    tmp_dat_dgc$Type <- "IM + DGC"

    ## Normal->IM->GC看基因的表达变化情况
    tmp_dat_gc <- tmp_dat
    tmp_dat_gc$Class <- ifelse( tmp_dat_gc$Class %in% c("IGC" , "DGC") , "GC" , tmp_dat_gc$Class )
    tmp_dat_gc$Class <- factor( tmp_dat_gc$Class , levels = c( "Normal" , "IM" , "GC" ) , order = T )
    tmp_dat_gc$Type <- "IM + GC"
    ## 一个既有IGC又有DGC，合并
    tmp_dat_gc <- tmp_dat_gc %>%
    group_by(Sample , Class) %>%
    summarize( Tpm = median(Tpm) , Type = unique(Type) )

    tmp_dat <- rbind( tmp_dat_igc , tmp_dat_dgc , tmp_dat_gc )
    tmp_dat$Type <- factor( tmp_dat$Type , levels = c("IM + GC" , "IM + IGC" , "IM + DGC") )

    out_name <- paste0( image_path , "/" , gene , ".TPM.pdf" )

    title <- paste0(
        "Gene : " , gene , "\n" 
    )

    plotTpm(tmp_dat = tmp_dat , out_name = out_name , title = title)

    out_name <- paste0( image_path , "/" , gene , ".TPM.tsv" )
    write.table(tmp_exp , out_name , row.names = F , sep = "\t" , quote = F)
}


##########################################################################################
## 分单个基因，看在一个人多个样本的改变情况
Hugo_Symbol <- gene
image_path <- out_path
dir.create(image_path , recursive = T)

describeGene(dat_expression = dat_expression , image_path = image_path )



